N <- 1000
XRandSd <- 1/sqrt(N)
XRandM <- 0
eta <- 0.5
MML <- round(6*N/eta)

b<-matrix(1,N,1)
bL<-sqrt(sum(b^2))
set.seed(300)

j1<-b
for(i in 1:(N*0.2)){
	j1[i]=-1
}

j2<-b
for(i in ((N*0.2)+1):(N*0.6)){
#for(i in 1:(N*0.2)){
	j2[i]=-1
}

R<-matrix(0,MML+2,2)#R(S1,T),R(S2,T)
q<-matrix(0,MML+2,1)#R(S1,S2)
L<-matrix(0,MML+2,2)#length of S1,S2. l*sqrt(N)

step <- function(x) ifelse(0 < x, 1, 0)

mlearn <- function(){
	for (i in 1:(MML+1)){
		L[i,1]<<-sqrt(sum(j1^2))
		L[i,2]<<-sqrt(sum(j2^2))
		R[i,1]<<-sum(j1*b)/(L[i,1]*bL)
		R[i,2]<<-sum(j2*b)/(L[i,2]*bL)
		q[i]<<-sum(j1*j2)/(L[i,1]*L[i,2])
		
		x<-matrix(rnorm(N,XRandM,XRandSd),N,1)
		uj1<<-sum(x*j1)
		uj2<<-sum(x*j2)
		j1 <<- j1 + eta*step(-uj2*uj1)*sign(uj2)*x
		j2 <<- j2 + eta*step(-uj1*uj2)*sign(uj1)*x
	}
	L[(MML+2),1]<<-sqrt(sum(j1^2))
	L[(MML+2),2]<<-sqrt(sum(j2^2))
	R[(MML+2),1]<<-sum(j1*b)/(L[(MML+2),1]*bL)
	R[(MML+2),2]<<-sum(j2*b)/(L[(MML+2),2]*bL)
	q[(MML+2)]<<-sum(j1*j2)/(L[(MML+2),1]*L[(MML+2),2])
}

system.time(#returns CPU Time Used
	#MAIN
	mlearn()
,FALSE)

setwd("~/Documents/lab2006/dat")	#for OSX
if(eta==0.1)write.table(cbind(L,R,q),format(Sys.time(),"PercML-%Y.%b.%d-%H%M%S-eta01.txt"), sep="\t",row.names=F,col.names=F)
if(eta==0.3)write.table(cbind(L,R,q),format(Sys.time(),"PercML-%Y.%b.%d-%H%M%S-eta03.txt"), sep="\t",row.names=F,col.names=F)
if(eta==0.5)write.table(cbind(L,R,q),format(Sys.time(),"PercML-%Y.%b.%d-%H%M%S-eta05.txt"), sep="\t",row.names=F,col.names=F)
if(eta==0.7)write.table(cbind(L,R,q),format(Sys.time(),"PercML-%Y.%b.%d-%H%M%S-eta07.txt"), sep="\t",row.names=F,col.names=F)
if(eta==0.9)write.table(cbind(L,R,q),format(Sys.time(),"PercML-%Y.%b.%d-%H%M%S-eta09.txt"), sep="\t",row.names=F,col.names=F)

#
plot(L[1:(MML+2),1]/sqrt(N),q[1:(MML+2),1],xlim=c(0.,1),ylim=c(-1.,1),type="l")#;par(new=T)
#
plot(L[1:(MML+2),1]/sqrt(N),R[1:(MML+2),1],xlim=c(0.,1),ylim=c(0.2,1),type="l")
par(new=T);plot(L[1:(MML+2),2]/sqrt(N),R[1:(MML+2),2],xlim=c(0.,1),ylim=c(0.2,1),type="l")
#
plot(L[1:(MML+2),1]/sqrt(N),R[1:(MML+2),1],type="l")
par(new=T);plot(L[1:(MML+2),2]/sqrt(N),R[1:(MML+2),2],type="l")
#
plot(1:(MML+2),R[,1],ylim=c(0,1),type="l")
par(new=T);plot(1:(MML+2),R[,2],ylim=c(0,1),type="l")
par(new=T);plot(1:(MML+2),q,ylim=c(0,1),type="l")
#
plot(1:(MML+2),R[1:(MML+2),1],ylim=c(0,1),type="l")
par(new=T);plot(1:(MML+2),R[1:(MML+2),2],ylim=c(0,1),type="l")
par(new=T);plot(1:(MML+2),q[1:(MML+2)],ylim=c(0,1),type="l")
